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We generalize the corner transfer matrix renormalization group, which consists of White's 
density matrix algorithm and Baxter's method of the corner transfer matrix, to three dimen- 
sional (3D) classical models. The renormalization group transformation is obtained through 
the diagonalization of density matrices for a cubic cluster. A trial application to 3D Ising model 
with m = 2 is shown as the simplest case. 
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§1. Introduction 

Variational estimation of the partition function has been one of the standard technic in statistical 
mechanics. For a two-dimensional (2D) classical lattice model defined by a transfer matrix T, the 
variational partition function per row is written as 

_ (V\T\V) 

x ~ (V\V) ' { } 

where \V) represents the trial state and (V\ is its conjugate; A is maximized when (V| and |V) 
coincide with the left and the right eigenvectors of T, respectively. In 1941 Klamers and WannieiB'i) 
investigated the Ising model, assuming that \V) is well approximated by a product of matrices 

V(...,i,j,k,l,...) = ...F^F^ k F M ... , (1.2) 

where i, j, k, I, etc., are the Ising variables, and i™ is a symmetric 2 by 2 matrix. The approximation 
is more accurate than both the mean-field and the Bethe approximations^ Baxter improved the 
trial state by introducing additional degrees of freedom E'Si) His variational state is defined as 

V(...,i,j,k,l,...)= £ ...F%F£F%..., (1.3) 
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where a,b,c,d, etc., are 2 n -state group spin variables. The tensor F % b contains 4 ■ 2 2n elements, 



and therefore it is not easy to optimize F*i — adjust the elements — so that A is maximized. 



ab 

_ - adjust the elements 
He solved the optimization problem by introducing the corner transfer matrix (CTM), and by 
solving self-consistent equations for the tensors.^ In 1985 Nightingale and Blote used Baxter's 
tensor product as a initial state in the projector Monte Carlo simulation for the Haldane systemi 
Baxter suggested an outline of generalizing his CTM method to 3D systems]!* however, the project 
has not been completed. 

Similar variational formulations have been applied to one-dimensional (ID) quantum systems, 
especially for S = 1 spin chains. The variational ground state |^} is given by a modified tensor 
product 

*(...,;,j,M,...)= E -4^4-, (1-4) 

...,a,b,c,d,e,... 

where the subscripts a, b, c, d, e, etc., are m-state group spin variables. Affleck, Lieb, Kennedy, and 
Tasaki (AKLT) showed that the ground-state of a bilinear-biquadratic S = 1 spin chain is exactly 
expressed by the tensor product with m = 2% The variational formulation has been generalized by 
Fannes et. al. for the arbitrary large m. §00 Now such ground state is called 'finitely correlated 
state© or 'matrix product state'© Quite recently Niggemann et. al© showed that the ground 
state of a 2D quantum systems can be exactly written in terms of a two-dimensional tensor product. 
Although 1^) in Eq.(1.3) does not look like \V) in Eq.(1.4), they are essentially the same. We can 
transform \V) into |\&) by obtaining A l ab from F % J b through a kind of duality transformation;!^ the 
opposite is also possible. 

The application of both |V) in Eq.(1.3) and in Eq.(1.4) are limited to translationally in- 
variant (or homogeneous) systems. In 1992 White established a more flexible numerical variational 
method from the view point of the real-space renormalization group (RG).00) Since his numerical 
algorithm is written in terms of the density matrix, the algorithm is called 'density matrix renor- 
malization group' (DMRG). White's variational state is written in a position dependent tensor 
productS 

*(..., i,j,k, 1,...)= £ ...A\ h B{ c C k cd D de ..., (1.5) 

...,a,b,c,d,e ; ... 

where A ab is not always equal to B ab , etc. This inhomogeneous property in |3>) makes DMRG 
possible to treat open boundary systems0 and random systems. Now the DMRG is widely 
used for both quantum!!' and classicalS'lHH) problems. Quite recently Dukelsky et. al. ana- 
lyzed the correspondence (and a small discrepancy) between DMRG and the variational formula 
in Eq.(1.4)@lll§3> 

The decomposition of the trial state into the tensor product tells us how to treat lattice models 
when we try to obtain the partition function. The essential point is to break-up the system 
into smaller pieces — like the local tensor F^ b in Eq.(1.3) or A ab in Eq.(1.4) — and reconstruct 
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the original system by multiplying them. According to this idea, the authors combine DMRG 
and Baxter's method of CTM, and proposed the corner transfer matrix renromalization group 
(CTMRG) method. 0© It has been shown that CTMRG is efficient for determinations of critical 
indices© or latent heats .0 

The purpose of this paper is to generalize the algorithm of CTMRG to three-dimensional (3D) 
classical systems. We focus on the RG algorithm rather than its practical use. In the next section, 
we briefly review the outline of CTMRG. The key point is that the RG transformation is obtained 
through the diagonalization of the density matrix. In §3 we define the density matrix for a 3D vertex 
model, and in §4 we explain the way to obtain the RG transformation. The numerical algorithm 
is shown in §5. A trial application with m = 2 is shown for the 3D Ising Model. Conclusions are 
summarized in §6. 

§2. Formulation in Two Dimension 



Fig. 1. Square cluster of a symmetric vertex model; the shown system is the example with linear dimension 2N = 6. 
The cross marks x show the boundary spins. 



The aim of CTMRG is to obtain variational partition functions of 2D classical models. Let us 
consider a square cluster of a 16-vertex model (Fig.l) as an example of 2D systems. We impose 
the fixed boundary condition, where the boundary spins shown by the cross marks point the same 
direction. In order to simplify the following discussion, we assign a symmetric Boltzmann weight 
Wijki = Wj kli = W ilk j to each vertex,© where i,j, k and / denote two-state spins (= Ising spins or 
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arrows) on the bonds. (Fig. 2(a)) 



Fig. 2. Boltzmann weight and transfer matrices. The dots represent spin variables inside the square cluster shown 
in Fig.l, and the cross marks represent the boundary spins, (a) Vertex weight W^ kl . (b) Half-row transfer matrix 
P a l- ( c ) Corner transfer matrix C ab . 

We employ two kinds of transfer matrices in order to express the partition function Z 2N of the 
square cluster with linear dimension 2N. One is the half-row transfer matrix (HRTM). Figure 
2(b) shows the HRTM P ab with length N = 3, where the subscripts a = (a 1 ,a 2 , . . . ,a N ) and 
b = b 2 , ■ ■ ■ , b N ), respectively, represent row spins — in-line spins — on the left and the right 
sides of the HRTM. We think of P a l as a matrix labeled by the superscript i. The other is the 
Baxter's corner transfer matrix (CTM),Bi'@ ) that represents Boltzmann weight of a quadrant of 
the square. Figure 2(c) shows the CTM C ab with linear dimension N = 3. The partition function 
Z 2N is then expressed as 

Z 2N = Trp = TrC\ (2.1) 

where p ab = (C 4 ) ab is the density matrix. From the symmetry of the vertex weight W^jy, the 
matrices P a l, C ab , and p ab are invariant under the permutation of subscripts. 

There are recursive relations between W, P, and C. We can increase the length of HRTM by 
joining a vertex 

^ = EV«i- (2-2) 
k 
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Fig. 3. Extensions of (a) the HRTM (Eq.(2.2)), and (b) the CTM. (Eq.(2.3)) 



where the extended row-spins are defined as a = (a, I) = (a 1} a 2 , ■ ■ ■ , a N , I) and b = (b,j) = 
(b 1 , b 2 , ■ ■ ■ , b N ,j). (Fig. 3(a)) In the same manner, the area of CTM can be extended by joining two 
HRTMs and a vertex to the CTM 

C ~al= E W tjkl P d iPjC cd , (2.3) 

cd kj 

where the extended row-spins a and b are defined as a = (a, I) = (a 1; a 2 , . . . , a N , I) and b = (6, i) = 
(&!, b 2 , ■ ■ ■ , a N ,i). (Fig. 3(b)) In this way, we can construct HRTM and CTM with arbitrary size N 
by repeating the extension Eqs.(2.2) and (2.3). 

It should be noted that the matrix dimension of both C ab and P a l increases very rapidly with their 
linear size N. The fact prevents us to store the matrix elements of C ab and P a l when we numerically 
calculate the partition function Z 2N - This difficulty can be overcomed by compressing CTM and 
HRTM into smaller matrices via the density matrix algorithmJlIHil) where the RG transformation 
is obtained through the diagonalization of the density matrix p ab = (C 4 ) ab . Let us consider the 
eigenvalue equation for the density matrix 

J2p ab A% = \ a A2, (2.4) 

b 

where X a is the eigenvalue in decreasing order Ai > A2 > • • • > 0, and A a = {Af, A 2 , . . .) T is the 
corresponding eigenvector that satisfies the orthogonal relation 

(A«,A' 3 )=^^ J 4f = ^. (2.5) 

a 
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It has been known that X a rapidly approaches to zero with respect to aJj'EJ) and that we can 
neglect tiny eigenvalues from the view point of numerical calculation. We consider only m numbers 
of dominant eigenvalues in the following; the greek indices run from 1 to m. The number m is 
determined so that Y^a=i becomes a good lower bound for the partition function Z 2N = Tr p. 
Equation (2.4) shows that for a sufficiently large m the density matrix p can be well approximated 

as 

m 

Pab ^J2 A a A b^- (2-6) 
a=l 

The above approximation shows that the m-dimensional diagonal matrix 

PaP = E A a A bPab = ^A/J (2.7) 
ab 

contains the relevant information of p; we can regard p as the renormalized density matrix. This 
is the heart of the density matrix algorithm: the RG transformation is defined by the matrix 
A = (A 1 , A 2 , . . . , A" 1 ) which is obtained through the diagonalization of the density matrix. 

As we have applied the RG transformation to the density matrix p, we can renormalize the CTM 
by applying the matrix A as 

C a p = Y, A a4C ab . (2.8) 

ab 

Since C ab and p ah have the common eigenvectors — remember that p = C 4 — the renormalized 
CTM is an m-dimensional diagonal matrix 

C = diag(o;i, u 2 , w m ) , (2.9) 

where u> a is the eigenvalue of the CTM that satisfies X a = oj^- In the same manner, we obtain the 
renormalized HRTM 

PJ p = J2A^ b P a i. (2.10) 

ab 

In this case is not diagonal with respect to a and /?; the RG transformation is not always 
diagonalization. 

We can extend the linear size of CTM and HRTM using Eqs.(2.2) and (2.3), and we can reduce 
their matrix dimension by the RG transformation in Eqs.(2.7) and (2.8). By repeating the extension 
and the renormalization, we can obtain the renormalized density matrix p and the approximate 
partition function Z 2N = Tt p for arbitrary system size N. This is the outline of the CTMRG. 

§3. Density matrix in Three Dimension 

In order to generalize the density matrix algorithm to 3D systems, we first construct the density 
matrix in three dimension. As an example of 3D systems, we consider a 64-vertex model, that is 
defined by a Boltzmann weight Wj,-^ TO „- (Fig.4(a)) In order to simplify the following discussion, 
we consider the case where W^ klrnn is invariant under the permutations of the two-state spins 
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i,j,k,l,m and n£$ As we have considered a square cluster in two-dimension, (Fig.l) we consider 
a cube with linear dimension 2N, where the boundary spins (on the surface of the cube) are fixed 
to the same direction. According to the variational formulation shown in §1, we first decompose 
the cube into several parts shown in Fig.4(b)-(d). 



Fig. 4. Parts of the cubic cluster with linear dimension 2JV: (a) Vertex weight Wj,-fc( mn . (b) The tensor P a l cd . (c) 
The tensor S? b Y . (d) Corner Tensor C XYZ . The cross marks x represent the boundary spins. 



The tensor P a l cd shown in Fig. 4(b) is a kind of three-dimensional HRTM. The superscript i 
represents the two-state spin at the top. The spin at the bottom is fixed, because it is at the 
boundary of the system. The subscript a represents the group of in-line spins a = (a l3 a 2 , ■ ■ ■ , a N ); 
b, c, and d are defined in the same way. From the symmetry of the vertex weight, P a l C( i is invariant 
under the permutations of subscripts. 

The tensor S^- shown in Fig. 4(c) does not have its 2D analogue; it is an array of vertices. 
The subscripts a and b represent in-line spins; other two sides are the boundary of the cube. The 
superscript X represents an by N array of spins on the square surface 



X 



11 



21 



12 



1 22 



J 1N 



\ 



\ X N1 



(3.1) 



~NN j 



where x NN is closest to the center of the cube, and Y is the spin array on the other surface; x ij 
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and are connected to the same vertex at the position {ij}- The tensor is invariant under the 
permutation of X and Y (S x - = S YX ), Du * ^ s n °t invariant for a and b {S XY / S XY ); S XY is 
equal to where Z = X T and W = Y T . 

Figure 4(d) shows the corner tensor C XYZ , which is a kind of three-dimensional CTM.i^ The 
superscripts are defined in the same way as Eq.(3.1). (The boundary spins on the surfaces of the 
original cube are fixed.) It should be noted that C XYZ is not equal to C where W = X T , 
because each surface has its own orientation. 



Fig. 5. Extensions of (a) P in Eq.(3.2), (b) S in Eq.(3.3), and (c) C in Eq.(3.5). 



Following the formulation in two-dimension, let us consider the size extension of P, S, and C. 
The length of P can be increased by joining a vertex (Fig. 5(a)) 

^ahcd = ^ijklmn ^abcd : (3-2) 
n 

where the extended in-line spins are defined as a = (a,j) = (a 1 , a 2 , ■ ■ ■ , a N , j), b = (b,k) = 
(6 1 ,6 2 > ■■■ ,b N ,k), c = (c,l) = (c 1} C2, . . .,c N ,l), and d = (d,m) = (d 1 ,d 2 , . . .,d N ,m). The lin- 
ear size of S can be increased by joining two P and a vertex (Fig. 5(b)) 

S^l = ^ijklmn ^abcd ^efgh ^ce^ (3-3) 

In ce 

where the extended in-line spins are defined as a = (a,j) = (a 1 ,a 2; • • • ,a, N , j), and b = (g,i) = 
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(<7i, g 2 , • • • , g N , i). The extended spin array X is denned as 

/ 

X = 



x ll 


X 1N 


fl 


X N1 


■ X NN 


In 


h ■ 


■ b N 


k 



(3.4) 



and Y is denned in the same way from the indices m, d, h and Y. The linear size of the corner 
tensor C can be increased by joining three P, three S, and a vertex (Fig.5(c)) 



C 



XYZ 



EV^ m p n p I p m nXT oYU oZV fiTUV 

/ j / j ijklmn r abcd ^efgh ■ r opqr'- > qd "ce '-'hr ^ 

Iran cd eh qr TUV 



(3.5) 



The extended superscripts X, Y, and Z are defined in the same way as Eq.(3.4). In equation (3.5) 
we have to take care of the orientation of the surfaces T, U, and V. 



Fig. 6. The density matrix Q in Eq.(3.8) is obtained by joining two corner tensors (Eq.(3.6)) to obtain the tensor 
D, and then joining four of them. (Eq. (3.8)) 



Now we can express the partition function E 2N of the cube with linear size 2N using the corner 
tensors. We first join two corner tensors (Fig.5) to obtain a symmetric matrix 

D (xu)(zv) = C XYZ C" YV , (3.6) 

Y 

where we regard the pair (ZV) as the column index of D, and (XU) as the row index. The tensor 
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C m is the mirror image of C: C^ YV = C uTYyT . The partition function E 2N is then expressed as 



~ 2N = TrD* = Y,Q (XU)(XU) , 
XV 

where the matrix Q is the forth power of D (Fig.6) 



(3.7) 



q(XU)(ZV) = D (XU)(AB) D (AB)(CD) D (CD)(EF) D (EF)(ZV) _ (33) 

(AB){CD)(EF) 

The matrix Q can be seen as a density matrix for the cube, because Tr Q is the partition function 
H 2Ar . By contracting two superscripts of Q, we obtain a density submatrix 



p xz = j2q(xu)(zu) 
u 



(3.9) 



which will be used for the RG transformation for the spin array. 

Let us consider a density submatrix p xz for the extended cube with size 2(N + 1), where X is 
the extended spin array Eq.(3.4); for a while we label the elements of Z as 



x n 


x 'lN 


f'l 




• • x 'nn 


I'n 


b\ . 


.. b' N 


k' 



(3.10) 



in order to define another density submatrix. By tracing out N by N + 1 variables of the extended 
density matrix p xz 



Pf9 



E 



p 



(3.11) 



b i= b 'i x lj = x 'v 



where / = (/ l5 . . . , f N , k) and g = (f' 1 ,...,f' N ,k'), we obtain another density submatrix for 
the extended in-line spins. In the same way, we obtain for the in-line spins of length TV - 
/ = (x 1N , . . . , x NN ) and g = (x' 1N , . . . , x' NN ) — by tracing out N — 1 by N variables of p xz in 
Eq.(3.9). 

§4. RG Algorithm in Three Dimension 

As we have done in §2, we obtain RG transformations by way of the diagonalizations of the 
density submatrices. We first consider the eigenvalue relation 

J2p XZ Ui = A*U$ , (4.1) 

z 

where we assume the decreasing order for A#. We keep first m! eigenvalues, (^ = 1, . . . , m') and 
neglect the rest of relatively small ones. We then obtain the RG transformation matrix U$ , that 
maps the spin array X to an m'-state block spin For example, the corner tensor C XYZ is 
renormalized as (Fig. 7(a)) 

C** e = ]T U X U Y U Z C XYZ . (4.2) 

XYZ 
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It should be noted that under the transpose of the spin array X — > X T the matrix U£ transforms 
as ±U£ according to the parity of the block spin <I>. 
Let us consider another eigenvalue relation 

Y,Pf 9 A t = X ' A 1 (4-3) 

9 

for the density submatrix p^ g , where / and g are in-line spins. We assume the decreasing order 
for A^ as before, and we keep m numbers of large eigenvalues, (ip = 1, . . . ,m) The matrix 
then represent the RG transformation for the in-line spin /. For example, P a \ cd is renormalized as 
(Fig.7(b)) 

Pjh6 = E PalcdK^AJA^ . (4.4) 

abed 

By using both U$ and Aj, we can renormalize S XY as (Fig. 7(c)) 

= E S$rj%4u$Ul . (4.5) 

abXY 

As a result of RG transformations, the tensors P abcd , S XY , and C XYZ are approximated as 

m 

Paled- E A a a A^ b A2A d P^ 5 (4.6) 

a/3jS=l 
m m' 

S% ~ E E ^Etf £tf S*J (4.7) 

a/3=l *$=l 
m' 

fl U x U Y U§C** e . (4.8) 
*$e=i 

For the models that have unique ground-state spin configuration, the above approximations become 
exact when T = and T = oo even for m = m! = 1. 



Now we can directly generalize the algorithm of CTMRG to 3D lattice models. The algorithm 
consists of the extensions for P a l cd (Eq.(3.2)), S x b Y (Eq.3.3)), and C XYZ (Eq.(3.5)), and the RG 
transformations Eqs.(4.2),(4.4) and (4.5). The procedure of the renormalization group is as follows: 

(1) Start from N = 1, where all the tensors can be expressed by the Boltzmann weight W- klmn : 

P abed = W iabedx' S aY = W aXbYxxi and C ^ ™ = W ZXYxxx' where the mark ' X ' represents 

the fixed boundary spin. 

(2) Join the tensors W i: j klmn , P a l cd , S XY , and C XYZ using Eqs.(3.2), (3.3), and (3.5), respectively, 
and obtain the extended ones P-i- d , S x ^ ', and C XYZ . (Increment N by one.) 

(3) Using the extended corner tensor C XYZ in Eq.(3.5), calculate the density matrix p xz via 
Eq.(3.9) and its submatrix pj_ in Eq.(3.11). 

(4) Obtain the RG transformation matrix U§ and A% using Eqs.(4.1) and (4.3), respectively. We 
keep m' states for Vf, and m states for a. 
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Fig. 7. The renormalized tensors (a) Pj 0yS in Eq.(4.4), (b) S*| in Eq.(4.5), and (c) <7** e in Eq.(4.2). The greek 
letters a, /3, 7, and S denote m-state renormalized in-line spins, and the capital ones <&, $, and O denote m'-state 
renormalized spin arrays. 



(5) Apply the RG transformations to the extended tensors to obtain P^p-yS (Eq.(4.4)), S 1 ** 
(Eq.(4.5)), and (Eq.(4.2)). 

(6) Goto the step (2) and repeat the procedures (2)- (5) for the renormalized tensors P, S, and C. 

Every time we extend the tensors in the step (2) the system size — the linear dimension of the 
cube — increases by 2. After the step (3) we can obtain the lower bound of the partition function 
by taking the trace of the density submatrix ^(N+i) = Trp xz = Trp^_. We stop the iteration 
when the partition function per vertex converges. Since the extended spin array X of the density 
matrix p xz contains the original (unrenormalized) spin variable, we can directly calculate the local 
energy and the order parameter© 

Let us apply the above algorithm to the 3D Ising model. The model is equivalent to the 64-vertex 
model whose vertex weight is given by 

W mmn = U ^U aj U ak U al U am U an , (4.9) 

a=±l 

where U ai is unity when a = i, and is e K + \Je 2K — 1 when a 7^ i. The parameter K denotes the 
inverse temperature J/k^T. For this model the initial conditions for step (1) are slightly modified 

as 

pi — TT TJ TT TJ TJ 

± abed ^ xi^ xo u x6 w xc u xd 
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Sat? — ^xX^xy^xa^xfe 

C XYZ = U xX U xY U xZ , (4.10) 

where 'x' represent the boundary Ising spin. (The modification is nothing but the change in 
normalizations.) We impose ferromagnetic boundary condition x = 1. As a trial calculation we 
keep only two states for both in-line spins (m = 2) and spin arrays (m' = 2); when m' = 2 the 
parity of the renormalized spin array \& in Eq.(4.1) is always even. Figure 8 shows the calculated 
spontaneous magnetization. Because of the smallness of m and m', the transition temperature T c 
is over estimated, where the feature is common to the Klamers-Wannier approximation .EP 



Fig. 8. Calculated spontaneous magnetization of the 3D Ising model when m = m' = 2. The arrow shows the true 



Compare to the CTMRG algorithm for 2D classical systems, the above RG algorithm for 3D 
system requires much more computational time. The reason is that after the step (2) the extended 
in-line spin / becomes 2m-state, and the extended spin array X becomes 2m 2 m'-state; in order to 
obtain p xz in the step (3) we have to create a matrix D^ xu ^ zv ' by Eq.(3.6), whose dimension is 
4m 4 m /2 . For the simplest (non-trivial) case m = m' = 2 the dimension is already 256. 

§5. Conclusion and discussion 

We have explained a way of generalizing the RG algorithm of CTMRGHii to 3D classical 
models, focusing on the construction of the density matrix from eight corner tensors. The RG 
transformations are obtained through the diagonalizations of the density matrices. As far as we 
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know it is the first generalization of the infinite-system density matrix algor ithmBH) to 3D classical 
systems. 

From the computational view point, the calculation in 3D is far more heavy than that of CTMRG 
in 2D; we have to improve the numerical algorithm in 3D for realistic use. What we have done is to 
approximate the eigenstate of a transfer matrix in 3D as a two-dimensional product of renormalized 
tensor P (Eq.(4.4)); the most important process is to improve the tensor elements in P so that 
the variational partition function is maximized. The improvement of tensor product state for ID 
quantum system proposed by Dukelsky et. al. JHHHt where their algorithm does not explicitly 
require the density matrix, may of use to reduce the numerical effort in three dimension. 

The authors would like to express their sincere thanks to Y. Akutsu, M. Kikuchi, for valuable 
discussions. T. N. is greatful to G. Sierra about the discussion on the matrix product state. K. O. 
is supported by JSPS Research Fellowships for Young Scientists. The present work is partially 
supported by a Grant-in- Aid from Ministry of Education, Science and Culture of Japan. Most of 
the numerical calculations were done by NEC SX-4 in computer center of Osaka University. 
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